## Copyright (C) 2001 Paulo Neis
##
## This program is free software; you can redistribute it and/or modify
## it under the terms of the GNU General Public License as published by
## the Free Software Foundation; either version 2 of the License, or
## (at your option) any later version.
##
## This program is distributed in the hope that it will be useful,
## but WITHOUT ANY WARRANTY; without even the implied warranty of
## MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
## GNU General Public License for more details.
##
## You should have received a copy of the GNU General Public License
## along with this program; If not, see <http://www.gnu.org/licenses/>.
## 

## usage: [Zz, Zp, Zg] = ncauer(Rp, Rs, n)
##
##Analog prototype for Cauer filter.
##[z, p, g]=ncauer(Rp, Rs, ws)
##Rp = Passband ripple
##Rs = Stopband ripple
##Ws = Desired order
##
## References: 
##
## - Serra, Celso Penteado, Teoria e Projeto de Filtros, Campinas: CARTGRAF, 
##   1983.
## - Lamar, Marcus Vinicius, Notas de aula da disciplina TE 456 - Circuitos 
##   Analogicos II, UFPR, 2001/2002.

## Author: Paulo Neis <p_neis@yahoo.com.br>

function [zer, pol, T0]=ncauer(Rp, Rs, n)

## Cutoff frequency = 1:
wp=1;

## Stop band edge ws:
ws=__ellip_ws(n, Rp, Rs);

k=wp/ws;
k1=sqrt(1-k^2);
q0=(1/2)*((1-sqrt(k1))/(1+sqrt(k1)));
q= q0 + 2*q0^5 + 15*q0^9 + 150*q0^13; %(....)
D=(10^(0.1*Rs)-1)/(10^(0.1*Rp)-1);

##Filter order maybe this, but not used now:
##n=ceil(log10(16*D)/log10(1/q))

l=(1/(2*n))*log((10^(0.05*Rp)+1)/(10^(0.05*Rp)-1));
sig01=0; sig02=0;
for m=0 : 30
	sig01=sig01+(-1)^m * q^(m*(m+1)) * sinh((2*m+1)*l); 
end
for m=1 : 30
	sig02=sig02+(-1)^m * q^(m^2) * cosh(2*m*l); 
end
sig0=abs((2*q^(1/4)*sig01)/(1+2*sig02));

w=sqrt((1+k*sig0^2)*(1+sig0^2/k));
#
if rem(n,2)
	r=(n-1)/2;
else
	r=n/2;
end
#
wi=zeros(1,r);
for ii=1 : r
	if rem(n,2)
		mu=ii;
	else
		mu=ii-1/2;
	end
	soma1=0;
	for m=0 : 30
		soma1 = soma1 + 2*q^(1/4) * ((-1)^m * q^(m*(m+1)) * sin(((2*m+1)*pi*mu)/n));
	end
	soma2=0;
	for m=1 : 30
		soma2 = soma2 + 2*((-1)^m * q^(m^2) * cos((2*m*pi*mu)/n));
	end
	wi(ii)=(soma1/(1+soma2));
end
#
Vi=sqrt((1-(k.*(wi.^2))).*(1-(wi.^2)/k));
A0i=1./(wi.^2);
sqrA0i=1./(wi);
B0i=((sig0.*Vi).^2 + (w.*wi).^2)./((1+sig0^2.*wi.^2).^2);
B1i=(2 * sig0.*Vi)./(1 + sig0^2 * wi.^2);

##Gain T0:
if rem(n,2)
	T0=sig0*prod(B0i./A0i)*sqrt(ws);
else
	T0=10^(-0.05*Rp)*prod(B0i./A0i);
end

##zeros:
zer=[i*sqrA0i, -i*sqrA0i];

##poles:
pol=[(-2*sig0*Vi+2*i*wi.*w)./(2*(1+sig0^2*wi.^2)), (-2*sig0*Vi-2*i*wi.*w)./(2*(1+sig0^2*wi.^2))];

##If n odd, there is a real pole  -sig0:
if rem(n,2)
        pol=[pol, -sig0];
end

##
pol=(sqrt(ws)).*pol;
zer=(sqrt(ws)).*zer;

endfunction
